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1.  Summary 

This  project  involves  the  investigation  of  the  wave  propagation  in  nonlinear  composite 
media.  Works  were  mainly  concentrated  on  obtaining  analytical  and  numerical  solutions  of  the 
(2+l)-D  nonlinear  Schrddinger  equation,  which  accurately  models  nonlinear  wave  propagation. 
This  effort  demonstrated  theoretically  that  we  can  form  a  composite  medium  from  a  nonlinear 
medium  with  optical  vortex  solitons  (OVS).  These  spatial  solitons  actually  induce  ideal 
waveguides  within  the  nonlinear  medium.  What  is  more,  the  characteristics  of  these  waveguides 
may  be  dynamically  controlled  with  all-optical  means.  Through  computer  simulations,  we 
showed  that  a  waveguide  modulator  whereby  an  optical  vortex  soliton  induced  a  cylindrical 
waveguide  within  a  self-defocusing  medium  and  another  beam  beam  was  used  to  modulate  the 
refractive  index  within  a  small  region  of  the  waveguide.  On  the  fundamental  understanding  of 
dark  solitons  with  two  transverse  degrees  of  freedom,  we  have  studied  polarized  OVS’s.  We 
have  found  that  only  circularly  polarized  OVS’s  are  stable  and  other  polarized  OVS’s  are 
susceptible  to  long  period  of  spatial  perturbation.  What  is  remarkable  here  is  that  an  unstable 
polarized  OVS  exhibits  striking  decay  dynamics.  To  perform  the  above  research,  we  have 
developed  a  program  for  solving  vector  nonlinear  Schrddinger  Equation.  The  results  of  this 
project  have  been  shared  with  Dr.  G.  A.  Swartzlander,  Jr.,  O.  N.  T.  Postdoctoral  Fellow,  and  Dr. 
A.  J.  Campillo,  Senior  Research  Scientist,  both  of  the  Naval  Research  Laboratory  and  have  been 
presented  in  various  national  conferences. 


2a.  Polarized  optical  vortex  solitons 

Linear  vortices  are  already  known  to  play  a  fundamental  role  in  the  scattering  of  radiation 
and  in  waveguides,  where  the  vortex  is  characterized  by  the  separable  function  of  the  azimuthal 
coordinate:  exp(iM(t)),  and  M  =  +\,+l,  •••  is  the  so-called  topological  charge.  Cylindrical 
symmetry  allows  various  polarization  modes,  e.  g.,  transverse-electric,  transverse-magnetic,  and 
circular.  One  may  expect  nonlinear  optical  vortices,  which  also  have  cylindrical  symmetry,  to 
exhibit  interesting  polarization  dynamics  owing  to  the  symmetry-breaking  effect  of  an 
anisotropic  light-induced  refractive  index  change.  We  have  discovered  that  out  of  six 
characteristically  different  polarized  vortex  states,  only  the  circular  polarized  waves  are  stable. 
What  is  also  remarkable  is  that  the  decay  dynamics  of  the  other  five  modes  exhibit  two  striking 
phenomena:  (1)  the  initial  vortex  core  vanishes  and  then  re-emerges  with  the  opposite 
topological  charge,  and  (2)  additional  pairs  of  vortices  are  generated  (while  preserving  the  net 
topological  charge). 

The  propagation  of  polarized  light  in  a  nonlinear  refractive  medium  is  described  by  the 
paraxial  wave  equation, 

— >  f  1  ~~^NL 

i2kodE/dz+  d^/dx^ +  d^/dy^  E  +  klD  /no  =  0  (la) 

where  no  is  the  linear  refractive  index  (^sumed  isotropic),  ko  is  the  wave  number  inside  the 
medium  directed  along  the  z-axis,  and  E  is  the  electric  field.  The  electric  displacement  for  a 
Kerr  nonlinear  medium  is  given  by, 

—^NL  — ^  ^ 

D  =2n2Zo[AiE-E  )E+B{E'E)E  ]  (lb) 

where  n2  is  the  nonlinear  refractive  index  coefficient  for  linearly  polarized  light  field,  and,  for  a 
lossless  medium,  B  =  1  -  A  is  a  measure  of  the  light-induced  anisotropy.  In  a  self-defocusing 
medium,  n2  is  negative,  an^^the  intensity-dependent  refractive  index  is  given  by  n=7^o“^^ 
where  An  =  \n2\  E(x,y,z)' E  (x,y,z). 

— ^ 

It  is  naturally  convenient  to  decomposed  E  into  two  circularly  polarized  components: 
E  =  Vix,y)exp{-i2Z),  V=  M^(x,y)e^ M^(x,y)  4  >  where  4  =  2“’^^  (f  - /y )  and 

ei  =  2~^^^  {x  +  iy)  are  the  unit  vector  for  the  right  {R)  and  left  (L)  circularly  polarized 
components  respectively,  E^  and  Ei  are  the  values  of  the  background  field  amplitudes  of  each 
circularly  polarized  component,  x  and  y  are  transverse  unit  vectors,  and  are  the  normalized 
amplitude  profile  of  the  field,  Z=z/L\\  is  the  normalized  propagation  distance 
[L||  =2ko*  («o/^^~)],  and  An^  is  the  index  change  attributed  to  the  background  field.  We 
ignore  small  field  components  in  the  £  direction  that  occurs  for  transverse  magnetic  waves. 
From  Eq.  (la)  and  (lb),  two  coupled  nonlinear  Schrodinger  (NLS)  equations  are  obtained: 

2ui^  +  iduii/dZ  +  V^U{i-2  (1-B) | m/j  | +  (1+B) | | (2a) 
2u[^  +  idui^/dZ  +  ViUi-2  (1-B) | | +  (1-t-B)  1  |  ui  =  0  (2b) 

where  V  i=d^/dX^  r\ii=Ef{/ \E^  \  and  r\i  =  Ei/  \E^  \  give  the  fraction  of  energy  in 

each  polarization  state  (1£„|=  {Ej^  +  Ej}^''^),  and  (X,y)  =  (x/L^,y/Ly)  are  the  normalized 
transverse  coordinates  (L|=2/:5*  (hq/ An^y^^). 

The  scalar  NLS  equation  may  admit  vortex  solutions  of  the  form  u^(p,(l))  =  \|/(p)exp(/Af({)), 
where  p  =  (X^  +  y^)'^^  and  ^  is  the  azimuthal  coordinate  in  the  transverse  plane.  Now  we  can 
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inquire  whether  polarized  waves  allow  stable  vortex  modes  in  a  self-defocusing  medium.  An 
arbitrary  vortex  field,  V,  may  be  composed  by  superimposing  a  pair  of  right-  and  left-handed 
polarized  vortices  of  topological  charges  M  and  N,  e.  g.,  V  =  From  this 

infinite  selection  of  modes,  we  examine  the  most  fundamental,  namely,  those  with  unit  charge 
(M,  N  =  +  \)  whose  polarization  components  are  either  equal  in  amplitude  {Er=Ei)  or 
dominated  by  one  (e.g,.  ^  =  0).  This  allows  us  to  identify  ^ix  polarized  vortex  modes  of  special 
s^nificance:  “linear”  “circular”  Fqv  or  “radial” 

Y^ad  =  =  2'/2£;jVi/(p)r,  “radial-compliment”  =  “<!>” 

=  Ef^iuiei^- u_iei)  =  2’^^exp(-/7t/2)£^\|/(p)(l),  and  “(|)-compliment”  = 

Eji{u^\ej^~  uiCi).  Note  that  for  circularly  polarized  waves  An^=An2Ef(i  whereas  other 
polarizations  under  investigation  have  En^  =  n2{E\+E\).  While  many  other  modes  involving 
higher  or  mixed  charges  and  elliptical  polarization  will  not  be  treated  here,  one  may  expect  that 
the  results  of  the  following  stability  analyses  can  provide  enough  insight  to  allow  at  least 
qualitative  descriptions  for  other  cases. 


Fig.  1.  Instability  growth  rates  as  a  function  of  normalized  modulation  period,  T.  The  calculated 
initial  growth  rate,  F,  is  shown  for  radial  perturbations  (dashed).  The  average  growth  rate,  y,  is 
calculated  numerically  (data  points)  and  fitted  (solid)  to  Eq.  (4)  with  a  single  parameter.  For 
both  cases  we  set  the  radial  vortex  core  size  to  tq  =  25|i7n,  the  modulation  amplitude  to  £  =  0.1, 
and  material  parameters  to  An ^/nQ  =  2y.\(r^  and  B  =  1/3.  The  initial  growth  rate,  F,  is  roughly 
2. 1  times  the  average  rate,  y. 


To  determine  whether  these  six  vortex  modes  are  stable,  we  perform  a  linear  stability 
analysis  for  each  of  them.  We  find  out  that  the  propagation  constant  F  for  these  perturbed  modes 
can  be  expressed  as  a  function  of  normalized  spatial  frequency  ki=2%/T  or  normalized  spatial 
period  T (in  physical  unit,  the  period  is  r  =  TLj^): 

T=+iki{ki+2AY^^,+i{2  +  k\Y^-  for  circular  polarization  (3a) 

and  F =+iki{A  +  k\Y^^ ,+ki{AB -k\Y^^  for  the  remaining  five  polarizations.  (3b) 

Note  that  only  the  last  expression  in  Eq.  (3b)  may  be  real  under  certain  conditions.  Therefore 
the  circular  polarization  cases  are  stable  to  radial  perturbations,  owing  to  the  fact  that  F  is  always 
imaginary.  This  is  not  a  surprise  since  the  NLS  equation  for  the  circular  polarization  case  is 
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identical  in  form  to  the  scalar  NLS  equation,  the  latter  has  already  been  found  to  be  stable  (both 
in  theory  and  practice). 

On  the  other  hand,  the  remaining  five  polarization  states  admit  an  exponential  growth  of  the 
perturbation.  This  may  have  been  expected  since  it  is  known  that  elliptically  polarized  plane 
waves  are  also  unstable  under  cross-phase  modulation.  From  the  last  expression  of  Eq.  (3b),  we 
notice  that  the  vortex  is  unstable  whenever  ki<kcr  =  2B  (i.  e.,  whenever  the  physical  period  of 
modulation  exceeds  t^r  =  7cL|5"*^^).  The  growth  rate  of  instability  in  the  last  expression  of  Eq. 
(3b)  can  be  written  in  terms  of  T  and  7^/- 

T=AB{\-Tcr/Tf'^Tcr/T.  (4) 

The  growth  of  the  instability  peaks  at  Topt  =  ^T,r  =  l.l,  as  depicted  in  Fig.  1,  is  a  plot  of  Eq.  (4) 
for  B  =  1/3.  Since  perturbations  used  in  the  linear  stability  analysis  preserve  the  cylindrical 
symmetry  of  the  initial  vortices  and  maintain  zero  field  at  the  center  of  the  vortices,  they  may  not 
be  the  dominant  mode  causing  the  strongest  instability. 

We  have  used  numerical  methods  to  perform  a  more  rigorous  nonlinear  instability  analysis 
which  allows  not  only  the  determination  of  growth  rates,  but  also  qualitative  information  on  the 
decay  process.  For  numerical  convenience,  we  have  examined  a  sinusoidal  perturbation  along 
the  x-axis: 


uyi/  +  esin(/:j^X) 


eR+EL 


u^M-EsmikiX) 


4 


(5) 


where  we  set  the  amplitude  of  modulation  to  8=0.1,  and  M  =J-1.  The  perturbation  modulates 
each  of  the  circular  components  with  the  same  frequency  ki,  but  n  out  of  phase. 

The  decay  dynamics  from  these  modulations  have  been  investigated  by  numerical  solving 
Eq.  (2)  with  the  beam  propagation  method.  A  vortex  (with  core  radius  /■o=Lj/V2  =25)i?n) 
having  one  of  the  five  polarizations  under  consideration  is  subjected  to  periodic  modulations 
(t  =  1.0mm)  on  a  512x512  transverse  grid,  and  a  moderate  nonlinearity  of  An^/nQ=2.Qx\0~^  is 
assumed  in  a  medium  with  B  =  1/3  at  a  wavelength  of  1  \Lm.  Under  these  conditions  the 
nonlinear  scaling  size,  Lj^,  corresponds  to  5  pixels  in  our  computations.  The  dynamic  features 
for  each  of  the  five  polarizations  are  very  similar.  Therefore,  we  present  here  only  calculations 
for  azimuthal  polarization,  V,),.  The  intensity  and  phase  evolution  are  shown  in  Figs.  2  and  3, 
respectively,  at  different  propagation  distances.  The  gray-scale  palette  maps  small  values  of 
intensity  or  phase  (modulus  2k)  into  dark  tones;  a  logarithmic  (linear)  scale  is  used  for  intensity 
(phase).  These  cross-sectional  images  span  only  12%  of  the  computation  grid  area.  The  left- 
handed  circular  component  resembles  the  corresponding  right  circular  component  after  it  is 
rotated  in  the  x-y  plane  by  90^^  about  the  center  of  the  initial  vortex;  thus  it  is  not  shown. 

As  expected,  these  figures  show  that  the  perturbation  has  clearly  destroyed  the  circular 
symmetry  of  the  vortex  located  at  the  origin  (x=y=0).  Remarkable  and  unexpected  dynamics 
occurs  as  the  beam  propagates.  For  example,  in  (B)  of  Figs.  2  and  3,  the  vortex  core  has 
completely  vanished!  Instead,  a  k  phase  front  resembling  a  dark  soliton  segment  is  formed  in 
the  vicinity  of  the  origin.  (This  may  be  interpreted  as  a  topological  vortex  that  osculates  the  x-y 
plane,  and  thus,  does  not  actually  violate  the  conservation  of  charge  principle.)  Also  note  that 
although  the  vortex  itself  has  disappeared,  the  global  phase  continues  to  have  the  same  sense  of 
rotation  (counter-clockwise).  As  the  beam  proceeds  further  to  (C)  in  Figs.  2  and  3,  three  vortices 
develop,  with  the  central  one  at  the  origin  having  the  opposite  topological  charge  of  the  initial 
vortex,  i.e.,  its  sense  of  rotation  has  reversed!  The  two  "new"  vortices  have  the  same  charge,  and 
thus,  the  net  vorticity  (EM,)  is  the  same  as  the  initial  condition.  (This  may  be  interpreted  as  a 
single  vortex  filament  or  "string"  passing  through  the  transverse  plane  three  times.)  The 
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complexity  appears  to  develop  more  quickly  as  more  vortices  are  generated,  as  evidenced  in 
frame  (D).  In  this  case,  the  vortex  filament  osculates  the  x-y  plane  in  several  locations;  a  short 
distance  further,  seven  vortices  appear.  As  the  beam  proceeds  further,  one  may  expect  that 
nonlinear  optical  turbulence  will  eventually  develop. 


Fig.  2.  Transverse  intensity  profiles  of  the  right-hand  circularly  polarized  component  of  a 
transverse  electric  vortex  mode  (core  radius  Ro-2S\X]n,  nonlinearity  =  2.0x10“^) 

perturbed  by  a  modulation  with  an  amplitude  10%  of  the  beam  intensity  and  a  period  of  1mm  at 
Z  =  0  (0  cm),  5.9  (9.4  cm),  7.5  (1 1.9  cm),  and  11  (18.2  cm)  in  (A),  (B),  (C)  and  (D)  respectively. 


Fig.  3.  Same  as  Fig  3,  except  phase  is  shown  (black  corresponds  to  zero  phase,  and  white  to  27t). 
The  global  vorticity  has  a  counter-clockwise  sense.  The  vortex  core  is  absent  in  (B)  and  (D); 
instead  7t  phase  fronts  exist  and  the  corresponding  intensity  profiles  show  dark  soliton  segments. 

In  general,  the  vector  vortex  dynamics  may  be  interpreted  as  two  co-linear  vortex  filaments 
of  opposite  circular  polarization  that  unravel  (owing  to  cross-phase  modulation)  as  they 
propagate  through  the  self-defocusing  medium.  If  one  of  the  polarized  filaments  is  absent  or  if 
the  nonlinear  coupling  is  "turned  off,  stable  propagation  is  expected  from  the  scalar  theory'.  If 
both  polarization  modes  are  present,  their  coupling  effect  will  depend  on  the  relative  power  in 
each  mode,  as  well  as  the  material  parameter,  B.  For  a  given  value  of  B,  we  expect  the  strongest 
coupling  when  the  two  modes  have  equal  strength,  i.  e.,  the  five  unstable  modes  described  above. 
Based  on  these  results,  one  may  expect  spontaneous  vortex  formations  that  appear  as  closed 
loops  in  three  dimensions,  which  may  occur  when  a  plane  wave  of  arbitrary  polarization  (though 
not  circular)  is  perturbed  at  different  periods  along  both  the  x-  and  y-axes. 

Possible  impact  of  these  research  results  is  a  potential  scheme  for  making  a  polarization 
sensitive  all-optical  switch.  For  example,  a  light-induced  waveguide  is  formed  in  a  right 
circularly  polarized  beam  with  an  OVS.  Another  low  intensity  probe  beam  with  different  color 
can  be  guided  along  this  light-induced  waveguide  along  the  dark  core  of  the  OVS.  If  another 
gate  probe  with  the  same  frequency  but  different  different  polarization  is  introduced  to  the 
guiding  beam,  the  polarized  OVS  may  become  unstable.  As  a  result,  the  guided  probe  beam  will 
be  switched  off  owing  to  the  disintegration  of  the  light-induced  wavelength. 
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2b.  All-optical  modulator 

Besides  the  all-optical  switch  making  use  of  polarization  effect,  another  idea  is  to  exploit 
the  property  that  the  size  of  the  light-induced  fiber  (“light  pipe”)  changes  with  the  intensity  of 
the  bright  region  of  the  guiding  beam.  Since  an  OVS  has  a  core  radius  ro  =  LyV^  where 
L|  =  2/co^  (nQ/Arioo)^^^,  n^-n2  \  E^\supl  and  is  the  background  field,  the  radius  of  an  OVS 
is  proportional  to  1/|£„|.  By  squeezing  and  relaxing  the  “light  pipe”,  we  can  encode 
information  into  the  flow  of  the  wave.  Therefore,  we  can  achieve  this  by  modulating  the 
intensity  of  the  bright  region.  We  can  easily  introduce  modulations  to  the  guiding  beam  by 
sending  a  control  beam  to  it.  Consequently,  whatever  signal  is  present  in  the  control  beam  will 
be  transported  to  the  guided  light.  Fig.  4  displays  the  relations  among  the  guiding  beam  (with  an 
OVS),  the  guided  probe  beam  and  the  gate  beam  as  well  as  the  increase  in  owing  to  the 
gate  beam. 


+  Z 


>  Z 


Fig.  4.  All-optical  modulator  based  on  the  intensity- 
dependent  core  size  of  an  OVS.  The  core  size  as  well  as 
the  probe  beam  guided  by  the  OVS  inside  a  guiding  beam 
is  modulated  by  the  gate  signal 


Guided  beam 
lied  here 


Guiding  beam  with  OVS 

Fig.  5.  Pinch-off  of  the  light-induced 
waveguide  (OVS).  Observe  that 
after  passing  through  the  gate  beam, 
the  guided  beam  (the  upper  picture) 
radiates  and  dissipates  while  the 
OVS  in  the  guided  beam  (the  lower 
picture)  returns  to  original  size. 


Owing  to  the  reduction  in  the  the  dimension  of  the  waveguide,  portion  of  the  guided  beam 
must  radiate  out.  If  the  change  in  is  abrupt,  a  substantial  percentage  of  the  guided  beam 
can  escape  from  the  waveguide.  Our  prediction  is  confirmed  by  results  of  numerical  simulations 
of  a  guiding  beam  with  ro  =  25|l/?i  and  nominal  A/r  oo  =  2x  1 0“^ .  Fig.  5  shows  the  topview  (along 
the  propagation  direction)  of  a  guided  beam  with  intensity  0.001  of  the  OVS’s  background 
intensity  and  a  guiding  beam  with  an  OVS  (“light  pipe”)  in  one  of  our  computer  simulations 
where  a  gate  beam  is  directed  to  the  guiding  beam  perpendicularly.  The  gate  beam  causes  A/;^ 
to  increase  to  6x10“*^  The  cross  section  of  “light  pipe”  shrinks  abruptly  for  the  portion  of 
guiding  beam  under  the  influence  of  the  gate  beam,  and  relaxes  to  original  size  when  the  effect 
of  the  gate  beam  diminishes.  The  shrinkage  of  the  “light  pipe”  blocks  off  a  portion  of  the 
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power  of  the  guided  beam.  As  a  result,  the  guided  beam  is  invisible  after  traversing  the  gate 
beam.  This  is  analogous  to  the  pinch-off  effect  happening  in  the  channel  (“light  pipe”)  of  a 
field  effect  transistor. 

We  have  also  examined  the  dynamics  of  the  guided  pulse  under  different  levels  of  gate 
beam  power.  Fig.  6  shows  the  the  simulation  results  of  applying  various  gate  signal  intensities 
around  the  longitudinal  distance  between  9.1Zq  and  9.2Zq.  The  different  gate  intensities  result 
in  An^  varying  from  10“^  to  8x10“'*.  As  the  gate  power  increases,  the  guided  power  within  the 
longitudinal  region  enclosed  by  the  gate  beam  jumps  up.  After  the  longitudinal  region, 
background  intensity  of  the  guiding  beam  returns  to  nominal  value  and  so  does  the  OVS’s  core 
radius.  This  change  in  the  waveguide  geometry  leads  to  the  leakage  of  the  guided  radiation  and 
the  drop  in  guided  radiation  within  a  circular  aperture  of  radius  =  2rQ.  Fig.  6  demonstrates  that 
the  transmission  coefficient  of  the  guided  beam,  which  measures  the  amount  of  radiation  passing 
through  the  circular  aperture,  reaches  a  steady  and  lower  value  at  a  distance  of  9.4Zo,  i.e  0.2Zo 
away  from  the  gate  beam  region.  This  steady  value  of  transmission  coefficient  after  the  gate 
beam  decreases  almost  linearly  as  the  change  in  nonlinear  refractive  index  increases. 


Fig.  6.  Transmission  coefficient  of  the  guided  beam  through  a  aperture  with  radius  =2rQ  versus 
the  propagation  distance.  The  gate  beam  locates  at  the  distance  between  9.1Zo  and  9.2Zo. 


2c.  Propagation  program  for  vector  nonlinear  Schrodinger  equation 

We  have  extended  a  program,  which  was  developed  to  simulate  3  dimensional  nonlinear 
wave  propagation,  to  solve  the  vector  nonlinear  Schrodinger  equation  in  Eq.  (2).  We  have 
ported  the  program  to  Cray  C-90  supercomputer.  As  a  result,  it  takes  less  than  3  minutes  of  cpu 
time  to  simulate  optical  wave  propagation  to  a  distance  of  three  diffraction  length.  The  results  of 
each  run  are  usually  stored  as  many  frames  of  pictures  with  typical  resolution  of  512  x  512  pixels 
in  HDF  (Hierarchical  Data  Format)  which  is  developed  by  the  National  Center  for 
Supercomputing  Applications.  This  image  storage  format  reduces  the  storage  space 
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substantially.  Various  programs  have  been  implemented  for  processing  the  resulting  images  and 
format  conversion. 

The  core  of  the  nonlinear  propagation  code  is  based  on  the  split-step  method  (or  beam 
propagation  method).  This  numerical  technique  was  originally  used  to  solve  the  scalar  nonlinear 
Schrbdinger  equation  at  a  given  transverse  plane  by  splitting  each  numerical  step  into  a  linear 
half  step  without  nonlinearity  and  a  nonlinear  half  step  without  diffraction.  We  modify  it  here  to 
solve  Eq.  (2)  which  consists  of  two  coupled  nonlinear  Schrodinger  equations.  The  nonlinear 
phase  correction  is  computed  in  real  space  whereas  the  linear  phase  contribution  from  diffraction 
is  computed  in  the  Fourier  transform  plane.  Therefore,  the  electric  field  is  transformed  into 
Fourier  components  during  the  linear  step  and  is  transformed  back  into  real  space  during  the 
nonlinear  step.  These  transformations  are  carried  out  by  the  fast  2D  FFT  (Fast  Fourier 
Transformation)  on  the  Cray  C-90.  Since  the  version  of  FFT'  on  the  Cray  supercomputer  are 
highly  optimized,  results  can  be  obtained  from  simulations  using  the  code  in  a  very  short  time. 

The  listing  of  the  new  nonlinear  propagation  program  can  be  found  in  Section  3.  This  code 
can  be  run  on  other  UNIX  workstations  with  the  IMSL  mathematical  library  and  the  HDF  library 
installed. 
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3.  Listing  of  the  nonlinear  propagation  program 

We  list  here  the  beam  propagation  program  (called  non2beam.F)  for  solving  vector 
nonlinear  Schrbdinger  equation  in  Fig.  7.  Non2beam.F  can  simulate  nonlinear  propagation  of 
the  circularly  polarized  components  of  a  beam  with  arbitrary  polarization  or  nonlinear  co¬ 
propagation  of  2  beams  with  different  frequencies.  This  program  requires  an  input  data  file 
which  is  listed  in  Fig.  8  (called  non2beam.dat)  which  contains  data  for  the  nonlinear 
propagation.  On  a  Cray  C-90  supercomputer  with  UNICOS  system,  the  program  is  compiled  by 
issuing  the  following  command: 

cf77  -o  run  non2beam.F  -Idf 

where  the  “-Idf”  flag  is  for  the  libraiy  of  HDF  subroutines  (at  the  Pittsburgh  Supercomputing 
Center,  this  library  is  put  in  /usr/locaMib/libdf.a  which  is  in  the  automatic  search  path  of  cf77. 
Cray  supercomputers  at  other  locations  may  put  it  under  other  directories.)  On  a  typical  UNIX 
workstation,  the  program  can  be  compiled  by  issuing  the  following  command: 

f77  -O  -o  non2beam  non2beam.F  -Idf  -limsllO 

where  the  flags  “-Idf”  and  “-limsllO”  are  for  the  libraries  of  HDF  subroutines  and  IMSL 
subroutines,  e.g.  f77  looking  for  libraries  with  names  libdf.a  and  libimsllO.a  under  its  regular 
search  path.  After  compilation,  the  program  can  be  executed  with  the  following  command: 

non2beam  <  non2beam.dat »  non2beam.dat 

This  command  appends  output  information  to  the  input  file  input.dat.  Other  output  files  are 
specified  in  non2beam.dat. 


Fig.  7.  Listing  of  the  nonlinear  propagation  program  non2beam.F 


program  non2beain 
c 

c  Program' for  solving  3D  nonlinear  Schrodinger  equation 
c  with  beam  propagation  method  for  two  beam 
c 

c  define  polar  for  decomposed  a  beam  into  circular  components 
c 

parameter  (mx=512  ,  md=inx+l ,  ray=512 ) 

#ifdef  Cray 

parameter  (ntable=100+2*  (mx+my)  ,  nwork=4*mx*my) 

#else 

#ifdef  sgi 

parameter  {ntbl=30+my+mx,ncp=2*mx) 

#else 

parameter  (ntbl=15+4*mx, ncp=2*mx) 

#endif 

#endif 

complex  u {md, my) ,phX (mx) , phY (my) , cpxi , prox , proy , dem, optnon 
complex  phlin (md, my) ,  temp,  tmpo(8),  u2(md,my),  temp2 
#ifdef  Cray 
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real  table (ntable) ,  work(nwork) 

#else 

#ifdef  sgi 

complex  wffl{ntbl) 
real  cpy(ncp) 

#else 

real  wffl(ntbl),  wff2(ntbl),  cpy(ncp) 

#endif 

complex  work (mx, my) 

#endif 

integer  iflag 

character  filnam*18, filnam3*18,ch*20,ansl*l,ans2*l, filnam4*18 
character  filnam2*18,  filnam23*18,  filnam24*18 
c 

c  Default  value  for  some  of  the  inputs 
c 

data  wlam,wid, eps, sclx,  scly,  zl  /l.e-6,  l.e-4,  l.e-3,1.,  l.,0./ 
data  isig, ipr , nfrm  /-I, 5,1/  filnam3  /'ft.d'/  filnam4  /'ifr.hdf'/ 
data  zinit  /O./  refnO/1./  filnam23  /'ft2.d'/  filnam24  /'ifr2.hdf'/ 
c 

c  Skip  the  first  line  in  the  data  file 
c 

2  format (ilO) 

3  format (e20 . 13 ) 
read (5 , * ) 
read(5,l)  filnam 
read (5,1)  filnam2 
read (5,1)  ansi 

if  (ansi  .eq.  'y' )  then 
read (5,1)  ch 

if  (ch  .ne.  ''  ")  read(ch,2)  nfrm 
read (5,1)  ch 

if  (ch  .ne.  "  '' )  read(ch,3)  zinit 
else 

read ( 5 , * ) 
read ( 5 , * ) 
endif 

read (5,5)  e2n 
5  format (70x, e20 . 13 ) 
iflag  =  0 

if  (abs(e2n)  .It.  l.e-30)  iflag=l 
read (5,1)  ch 

if  (ch  .ne.  "  ")  read(ch,3)  wlam 
read (5,1)  ch 

■  if  (ch  .ne.  ''  ")  read(ch,3)  wid 
c 

c  Get  the  scale  factor  for  the  x  and  y  range 
c 

read (5,1)  ch 

if  (ch  .ne.  "  ")  read(ch,3)  sclx 
read (5,1)  ch 

if  (ch  .ne.  ''  ")  read(ch,3)  scly 
read (5,1)  ch 
1  format (70x, a) 

if  (ch  .ne.  "  ")  filnam3=ch 
#ifdef  Cray 

open(unit=3,  file=f  ilnam3  ,  status= 'unknown'  ,  access= '  direct '  , 
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1  recl=128,  forin='unformatted' ) 

#else 

open (unit=3 , f ile=f ilnam3 , status= 'unknown' , access= ' direct ' 
1  recl=16, form=' unformatted' ) 

#endif 

read (5,1)  ch 

if  (ch  .ne.  "  ")  filnam23=ch 
#ifdef  Cray 

open (unit=7 , f ile=f ilnam23 , status= 'unknown' , access= '  direct 
1  recl=128 , form= 'unformatted' ) 

#else 

open (unit=7 , f ile=f ilnam23 , status= ' unknown ' , access= ' direct 
1  recl=16 , form= 'unformatted' ) 

#endif 

pi  =  4 . *atan ( 1 . ) 
pi2  =  2 . *pi 
c 

c  e2n  --  n2*|E0P2/n0 
c  refnO  --  linear  refractive  index 
c  isig  --  sign  of  nonlinearity 

c  xrag  --  positive  range  of  x  in  meter,  [-xrag,  xrag] 
c  yrag  --  positive  range  of  y  in  meter,  [-yrag,  yrag] 
c  zrag  --  final  propagation  distance  z/zo 
c  wlam  --  wavelength  in  meter 
c  wid  --  beam  width  in  meter 
c  eps  --  error  control 

c  pratio  --  power  ratio  of  beam  2  to  beaml 

c  couple  --  coupling  coef.  between  beaml  and  beam2  (l+b)/a 
c 

read (5,1)  ch 

if  (ch  .ne.  "  ")  read(ch,3)  refnO 
read (5, 1)  ch 

if  (ch  .ne.  "  ")  read(ch,2)  isig 
read ( 5 , 1 ) ch 

if  (ch  .ne.  "  ")  read(ch,3)  zrag 
#ifdef  Cray 

open (uni t=2 , f ile=f ilnam, status= ' old' , f orm= ' unformatted '  , 

1  access= ' direct ', recl=128 ) 

#else 

open (uni t= 2 , file=f ilnam, status=  'old'  ,  form= 'unformatted' , 

1  access= ' direct ', recl=16) 

#endif 

read ( 2 , rec=l )  nx,  ny 

xrag  =  sclx*wid*0 . 5*sqrt (nx*pi) 

yrag  =  scly*wid*0 . 5*sqrt (ny*pi) 

wnum  =  2.*pi/  wlam 

xlen  =  xrag/wid 

ylen  =  yrag/wid 

zo  =  refn0*pi*wid**2/wlam 

zdist  =  zrag  *  0.25 

write (6, 6)  xrag, xrag,  yrag,  yrag, zdist, zo, wnum 
6  f ormat ( "  x/wid  range  [-",el4.7,"  m,  ”,el4.7,"  m]  ;  ", 

1  "y/wid  range  [-",el4.7,"  m,  '',el4.7,"  m]  "  /  “  z/4/Zo  is  " 
1  ,fl0.3,"  (Zo  '',el4.7,"  m)  ;  wavelength  number  -  ",el4.7 

1  ,"  1/m") 

read (5,1)  ch 

if  (ch  .ne.  "  ")  read(ch,3)  eps 
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ntot  =  nx*ny 
nx2  =  nx  +  2 

ny2  =  ny  +  2 

fnx  =  nx 

fny  =  ny 

nxd2  =  nx  /  2 
nyd2  =  ny  /  2 
cpxi  =  (0 . , 1 . ) 
delx  =  2.*xlen/nx 
dely  =  2 . *ylen/ny 
c 

c  delz  is  small  enough  so  that  the  linear  phase  <=  pi/2 
c 

delz  =  min ( 1 ./( float (nx) /sclx**2+float (ny) /scly**2 ), eps) 
iter  =  int (zdist/delz)  +  1 
write(6,20)  iter 

20  f ormat ( '  number  of  iterations  estimated  =  i5) 

read (5,1)  ans2 
if  (ans2  .eq.  'y')  then 
read (5,1)  ch 

if  (ch  .ne.  "  ")  filnam4  =  ch 
read (5,1)  ch 

if  (ch  .ne.  ''  ")  filnam24  =  ch 
c 

c  define  reference  level 
c 

rlv  =  0.5 
else 

read (5, *) 
endif 

read (5,1)  ch 

if  (ch  .ne.  "  " )  read(ch,2)  ipr 
majlp  =  iter/ipr 
if  (mod (iter, ipr)  .ne.  0)  then 
ma j Ip  =  ma j Ip  +  1 
iter  =  majlp*ipr 
endif 
c 

c  Get  beginning  distance  and  ending  distance 
c 

read (5,1)  ch 

if  (ch  .ne.  "  ")  read(ch,3)  zl 
read (5, 5)  pratio 
read (5, 5)  couple 

write (6, 27)  zl,  e2n,  pratio,  couple 
27  f  ormat  ( ''  zl  =  ",el4.7,  ”  nonlinearity  =  ",el4.7/ 

1  "  power2/powerl  =  ",el4.7,"  coupling  coef.  =  ",el4.7) 

zst  =  zl  *  0.25 
delz  =  zdist/iter 

write ( 6 , 25 )  nx, ny , iter , delz , ipr , eps 
25  f ormat ( lx, i5 , '  x  ',i5,'  x-y  dimension,  ', 

1  '  number  of  iteration  ',i5//'  dz  ',el4.7'  write  ', 

1  i5,'  frames  with  tolerance  <  ',el4.7) 

c 

c  compute  linear  propagation  in  x  direction 
c 

if  (if lag  .eq.  1)  then 
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prox  =  cpxi*delz*f loat (majlp) * (pi/xlen) **2 
else 

prox  =  cpxi*delz* (pi/xlen) **2*0 . 5 
endif 

phX(l)  =  1. 

phX(nxd2+l)  =  exp(prox* (nxd2) **2) 
do  3  0  i  =  2 ,  nxd2 

phX(i)  =  exp(prox* (i-1) **2) 

30  phX(nx2-i)  =  phX(i) 
c 

c  compute  linear  propagation  in  y  direction 
c 

if  (if lag  .eq.  1)  then 

proy  =  cpxi*float (majlp) *delz* (pi/ylen) **2 
else 

proy  =  cpxi*delz* (pi/ylen) **2*0 . 5 
endif 

phY(l)  =  1. 

phY(nyd2+l)  =  exp(proy* (nyd2) **2) 
do  40  i  =  2,  nyd2 
idx  =  i  -  1 

phY(i)  =  exp (proy* idx* *2) 

4  0  phY ( ny 2 - i )  =  phY ( i ) 
do  50  j  =  1,  ny 
do  50  i  =  1,  nx 

50  phlin(i,j)  =  phX(i) *phY( j ) 
zprop  =  0 . 

dntot  =1.  /  float (ntot) 
c 

c  recompute  phX  phY  for  later  use 
c 

prox  =  cpxi* (pi/xlen) **2 
phX(l)  =  0. 

phX(nxd2+l)  =  prox* (nxd2) **2 
do  53  i  =  2,  nxd2 

phX(i)  =  prox* (i-1) **2 
53  phX(nx2-i)  =  phX(i) 

proy  =  cpxi* (pi/ylen) * *2 
phY(l)  =  0. 

phY(nyd2+l)  =  proy* (nyd2 ) **2 
do  51  i  =  2,  nyd2 
idx  =  i  -  1 
phY(i)  =  proy* idx* *2 

51  phY(ny2-i)  =  phY(i) 
c 

c  get  initial  field  distribution 
c 

idx2  =  ntot* (nfrm-1) /8  +  1 
do  10  i  =  1,  ntot-7,  8 
idx2  =  idx2  +  1 

10  read  (2  ,  rec=idx2 )  (u  (mod  ( i+1:,  nx) +1 ,  ( i+k) /nx+1 )  ,  k=-l ,  6  ) 
close  (2) 

#ifdef  Cray 

open (uni t= 4 , f ile=f ilnam2 , status= ' old' , f orm= ' unformatted ' 
1  access= ' direct recl=128) 

#else 

open (unit =4 , f ile=f ilnam2 , status= 'old' , form= 'unformatted ' 
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1  access= ' direct recl=16) 

#endif 

idx22  =  ntot* (nfrm-1) /8  +  1 
do  310  i  =  1,  ntot-7,  8 
idx22  =  idx22  +  1 

310  read (4 , rec=idx22 )  (u2  (mod (i+k, nx) +1,  ( i+k) /nx+1) , k=-l , 6 ) 
close  (4) 

#ifndef  cray 
#ifdef  sgi 

call  cfft2di (nx,ny,wf fl) 

#endif 

#endif 

if  (ansi  .eq.  'n')  then 
c 

c  split2  routine  ’ 
c 

do  210  i  =  1,  ny 
do  215  3  =  1,  nxd2 
temp  =  u ( j , i ) 
temp2  =  u2 ( j , i) 
jnx2  =  j+nxd2 
u  ( j  ,  i)  =  u ( jnx2 , i) 
u2(j,i)  =  u2(jnx2,i) 
u2(jnx2,i)  =  temp2 
215  u(jnx2,i)  =  temp 
210  continue 

do  220  j  =  1 ,  nyd2 
do  225  i  =  1 ,  nx 
temp  =  u(i,j) 
temp2  =  u2 ( i , 3 ) 
jny2  =  3+nyd2 
u{i, j)  =  u(i, jny2) 
u2  (i , 3 )  =  u2 (i, jny2) 
u2(i,3ny2)  =  temp2 
225  u(i,3ny2)  =  temp 

220  continue 
c 

c  end  split2  routne 
c 

#ifdef  Cray 

call  cfft2d(l,nx,ny,l.,u, 1 ,md,u, 1 ,md, table, ntable, 

1  work,nwork) 

call  cf f t2d ( 1 , nx,ny, 1 . , u2 , 1 ,md, u2 , 1, md, table , ntable, 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cfft2d(-l,nx,ny,u,md,wffl) 
call  cf f t2d ( -1 , nx, ny , u2 , md, wf f 1) 

#else 

call  f  2t2d(nx,  ny  ,u,md,u,md,  wf  fl ,  wf  f  2  ,  work,  cpy) 
call  f 2t2d (nx, ny,u2  ,md,u2 ,md, wf f 1 , wf f2 , work, cpy) 

#endif 

#endif 

end  if 
c 

c  Propage  zl  distance  in  front  of  the  medium 
c 
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if  (abs(zst)  .gt.  l.e-30)  then 
do  70  i  =  1,  ny 
do  70  j  =  1 ,  nx 

u2(j,i)  =  u2(j,i)  *  exp{ {phX(j)+phY(i) ) *zst) 

70  u(j,i)  =  u(j,i)  *  exp ( (phX( j ) +phY{i) ) *zst) 
endif 

#ifdef  polar 

optconst  =  -2 . * (wnum*wid) **2*isig*e2n*refn0/ (1 . +couple) 

#esle 

optconst  =  -4 . * (wnuin*wid) **2*isig*e2n*refn0/ (1 . +couple) 

#endif 

phzinit  =  optconst*0 . 25*zinit 
c 

c  write  the  initial  fields 
c 

write (3 , rec=l)  nx,  ny 
write (7 , rec=l)  nx,  ny 
koo3  =  1 

if  (ans2  .eq.  'n')  then 
do  65  i  =  1,  ntot-7,  8 
koo3  =  koo3  +  1 

write  (7  ,  rec=koo3 )  (u2  (inod(i+k,nx)  +1 ,  ( i+k)  /nx+1 )  ,  k=-l ,  6 ) 

65  write (3 , rec=koo3 )  (u{mod(i+k,nx) +1, (i+k) /nx+1) , k=-l , 6) 
else 

#ifdef  Cray 

call  cfft2d(-l,nx,ny,dntot,u,  l,ind,u,  l,ind,  table, ntable, 

1  work,nwork) 

call  cfft2d(-l,nx,ny,dntot,u2, l,md,u2, l,md, table , ntable , 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cfft2d(l,nx,ny,u,ind,wffl) 
call  cscal2d(nx,ny,dntot,u,ind) 
call  cfft2d(l,nx,ny,u2,ind,wffl) 
call  cscal2d(nx,ny,dntot,u2,ind) 

#else 

call  f 2t2b (nx, ny , u,md, u, md, wf f 1 , wf f 2 , work, cpy ) 
call  f2t2b(nx,ny,u2,ind,u2,md,wffl,wf f2,work,cpy) 
do  279  i  =  1,  ny 
do  279  j  =  1,  nx 
u2(j,i)  =  dntot*u2 ( j , i) 

279  u(j,i)  =  dntot*u(j,i) 

#endif 

#endif 

call  hdf_iing2  ( filnam4  ,  phzinit ,  rlv,  u,  rnd,  nx,  ny,  0,work,ntot,  2) 
call  hdf_img2 ( f ilnain24 , phzinit , rlv, u2 , md, nx, ny , 0 , work, ntot , 2 ) 
#ifdef  Cray 

call  cf f t2d (1, nx, ny, 1 . ,u, 1 ,md,u, 1 ,md, table, ntable, 

1  work,nwork) 

call  cfft2d(l, nx, ny , 1 . , u2 , 1 ,md, u2 , 1 , md, table , ntable , 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cfft2d{-l,nx,ny,u,md,wffl) 
call  cfft2d(-l,nx,ny,u2,md,wffl) 

#else 

call  f2t2d (nx, ny , u, md, u, md, wf f 1 , wf f 2 , work, cpy) 
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call  f 2t2d (nx, ny, u2 ,md,u2,md,wffl,wf f2,work, cpy) 

#endif 

#endif 

endif 
energ  =  0 . 
energ2  =  0 . 
do  60  i  =  1,  ny 
do  60  j  =  1,  nx 

energ2  =  energ2  +  real (u2 ( j , i) *conjg (u2 ( j , i) ) ) 

.60  energ  =  energ  +  real (u { j , i) *conjg (u ( j , i) ) ) 
write (6,*)  zprop,  energ*dntot,  energ2*dntot 
c 

c  Start  the  split  step  method:  linear  part*nonlinear  part*linear  part 
c 

roptnon  =  optconst*delz 
optnon  =  cpxi*roptnon 
roptnon  =  float (majlp)  *  roptnon 
do  120  i  =  1,  ipr 
if  (iflag  .eg.  1)  then 
do  113  k  =  1 ,  ny 
do  113  jj  =  1,  nx 
u2(jj,k)  =  u2(jj,k)  *  phlin(jj,k) 

113  u(jj,k)  =  u(jj,k)  *  phlin(jj,k) 
zprop  =  zprop  +  delz*majlp 
else 

do  110  j  =  1,  majlp 
zprop  =  zprop  +  delz 
c 

c  linear  half  step 
c 

do  103  k  =  1 ,  ny 
do  103  jj  =  1,  nx 
u2(jj,k)  =  u2(jj,k)  *  phlin(jj,k) 

103  u(jj,k)  =  u(jj,k)  *  phlin(jj,k) 

#ifdef  Cray 

call  cfft2d(-l,nx,ny, dntot , u, 1 ,md, u, 1 ,md, table , ntable , 

1  work,nwork) 

call  cfft2d(-l,nx,ny, dntot , u2 , l,md, u2 , 1 ,md, table, ntable, 

1  work,nwork) 

#else 

ttifdef  sgi 

call  cf ft2d(l,nx,ny,u,md,wffl) 
call  cscal2d(nx,ny, dntot, u,md) 
call  cfft2d(l,nx,ny,u2,md,wffl) 
call  cscal2d(nx,ny, dntot, u2,md) 

#else 

call  f2t2b(nx,ny,u,md,u,md,wffl,wff2 , work, cpy) 
call  f 2t2b (nx,ny,u2 ,md,u2 ,md, wf f 1, wf f 2 , work, cpy) 
do  179  ii  =  1,  ny 
do  179  ji  =  1,  nx 
u2(ji,ii)  =  dntot*u2 ( ji, ii) 

179  u(ji,ii)  =  dntot*u( ji, ii) 

# endif 
# endif 
c 

c  nonlinear  full  step 
c 
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100 


do  100  k  =  1,  ny 
do  100  jj  =  1,  nx 
temp  =  u( j j ,k) *conjg(u( j j ,k) ) 
temp2  =  u2 { j j , k) *conjg(u2 ( j j , k) ) *pratio 
u2{jj,k)  =  u2 (j j ,k) *cexp(optnon* (temp2+couple*temp) ) 
100  u(jj,k)  =  u( j j ,k) *cexp{optnon* {temp+couple*temp2) ) 

c  dem  =  (l.jO.dO)  +  optnon*u{k) *conjg (u (k) ) 

c  100  u(k)  =  u(k)  *  ((-l.,0.)  +  (2.,0.)/dem)**2 
c 

c  linear  half  step 
c 

#ifdef  Cray 

call  cf f t2d (1 , nx, ny , 1 . ,u, 1 ,md,u, 1  ,md,  table , ntable , 

1  work.nwork) 

call  cf f t2d( 1 , nx, ny, 1 . , u2 , l,md, u2 , 1, md, table, ntable , 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cfft2d(-l,nx,ny,u,md,wffl) 
call  cfft2d{-l,nx,ny,u2,md,wffl) 


#else 


#endif 
#endif . 


call  f 2t2d (nx, ny , u, md, u, md, wf f 1 , wf f 2 , work, cpy ) 
call  f 2t2d (nx, ny, u2 ,md, u2 ,md, wf f 1 , wf f 2 , work, cpy) 


do  105  k  =  1,  ny 
do  105  jj  =  1,  nx 
u2(jj,k)  =  u2(jj,k)  *  phlin(jj,k) 
105  u(jj,k)  =  u(jj,k)  *  phlin(jj,k) 

110  continue 
endif 


c  output  field  distributions 


energ  =  0 . 
do  116  k  =  1 ,  ny 
do  116  jj  =  1,  nx 

116  energ  =  energ  +  real (u ( j j , k) *conjg (u ( j j , k) ) ) 

energ2  =  0 . 
do  316  k  =  1,  ny 
do  316  jj  =  1,  nx 

316  energ2  =  energ2  +  real (u2 ( j j , k) *conjg (u2 ( j j , k) ) ) 
if  (ans2  .eq.  'n')  then 
do  117  k  =  1,  ntot-7,  8 
koo3  =  koo3  +  1 

write (7 , rec=koo3 )  (u2 (mod(k+k0 , nx) +1 ,  (k+kO ) /nx+1 ) , k0=-l , 6) 

117  write (3 , rec=koo3 )  (u (mod (k+kO , nx) +1 , (k+kO) /nx+1 ) , k0=-l , 6 ) 
else 

#ifdef  Cray 

call  cfft2d(-l,nx,ny,dntot,u, l,md,u, l,md, table , ntable , 

1  work,nwork) 

call  cfft2d(-l,nx,ny,dntot,u2, 1 , md, u2 , 1 , md, table , ntable , 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cf ft2d(l,nx,ny,u,md,wf fl) 
call  cscal2d(nx,ny,dntot,u,md) 
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call  cfft2d(l,nx,ny,u2,md,wffl) 
call  cscal2d(nx,ny,dntot,u2,ind) 

#else 

call  f 2t2b (nx, ny,u,md,u,md, wf f 1 , wf f 2 , work, cpy) 
call  f2t2b(nx,ny,u2,md,u2,md,wffl,wf f2,work,cpy) 
do  379  ii  =  1,  ny 
do  379  ji  =  1,  nx 
u2(ji,ii)  =  dntot*u2 ( ji, ii) 

379  u(ji,ii)  =  dntot*u ( j i , ii) 

#endif 

#endif 

c 

c  Substract  a  reference  phase 
c 

refph  =  mod(phzinit+roptnon*float (i) ,pi2) 

refph2  =  mod (phzinit+roptnon*pratio* *  float ( i ), pi2 ) 

write (6, 129)  refph,  refph2 

129  f ormat ( '  max.  nonlinear  reference  phase  =  ',el4.7) 

call  hdf_img2 ( f ilnami , refph, rlv, u,md, nx, ny, i , work, ntot , 2 ) 
call  hdf_img2 (filnam24,refph2,rlv,u2,md,nx,ny, i, work, ntot, 2) 
#ifdef  Cray 

call  cfft2d(l,nx,ny,l.,u, 1 ,md,u, 1 ,md, table, ntable , 

1  work,nwork) 

call  cf f t2d (1, nx, ny, 1 . , u2 , 1 ,md, u2 , 1 , md, table, ntable, 

1  work,nwork) 

#else 

#ifdef  sgi 

call  cfft2d{-l,nx,ny,u,md,wffl) 
call  cf ft2d(-l,nx,ny,u2,md,wffl) 

#else 

call  f2t2d(nx,ny,u,md,u,md,wf fl,wf f2 , work, cpy) 
call  f2t2d (nx, ny , u2 ,md,u2 ,md,wf f 1 , wf f 2 , work, cpy) 

#endif 

#endif 

endif 

120  write (6,*)  4.*zprop,  energ*dntot,  energ2*dntot 
if  (ans2  .eq.  'y' )  then 
do  217  k  =  1,  ntot-7,  8 
koo3  =  koo3  +  1 

write (7 , rec=koo3 )  (u2 (mod(k+k0,nx) +1 , (k+kO) /nx+1) , k0=-l , 6 ) 

217  write (3 , rec=koo3 )  (u (mod (k+kO , nx) +1 , (k+kO ) /nx+1 ) , k0=-l , 6 ) 
endif 
c 

c  Propage  zdist+z2-zprop  distance  behind  the  medium 
c 

close  (3) 
close  (7) 
stop 
end 


Fig.  8.  An  example  of  the  input  file  non2beam.dat  in  which  the  data  from  the  user  are  bold¬ 
faced. 

DATA  FOR  NON2BEAM.F  enter  data  after  this  vertical  line  -->  | 

Name  of  the  file  containing  E-field  data  for  beaml  Igaus.d 

Name  of  the  file  containing  E-field  data  for  beam2  Igaus.d 


19 


n 


Is  the  above  file  produced  by  non2beam.f?  (y/n) 

Frame  #  to  use  if  the  previous  answer  is  y  (default  1) 

Initial  propagation  distance  if  the  previous  answer  is  y  (default  0) 
delta  n  (n2*E0**2/n0) 

Wavelength  in  meter  (default  l.e-6) 

Radial  size  wO  in  meter  (default  l.e-4) 

Resolution  for  x  range  (<1  for  high  resolution,  default  1) 
Resolution  for  y  range  (<1  for  high  resolution,  default  1) 

Name  of  the  file  containing  FT (field)  for  beaml  (default  ft.d) 

Name  of  the  file  containing  FT (field)  for  beam2  (default  ft2.d) 
Linear  refractive  index  nO  at  the  carrier  freq  (default  1) 

Sign  of  nonlinearity  (-1  or  +1)  (default  -1) 

Max.  Length  of  medium  in  z0==pi*n0*w0**2/  wavelength 
relative  error  (<<1)  (default  l.e-3) 

Output  images  in  hdf  format?  (y/n) 

If  y,  give  name  of  the  hdf  file  for  beaml  (default  ifr.hdf) 

If  y,  give  name  of  the  hdf  file  for  beam2  (default  ifr2.hdf) 

No.  of  frames  ipr  write  to  file  (default  5) 

Distance  of  the  beam  waist  in  front  of  the  medium  (in  ZO) 

Power  of  beaml  /  Power  of  beam2 
Coupling  coef.  between  the  two  beams 


|1. 

I 

|y 

I 

I 

1 

|0. 

|i. 

10. 


The  explanations  for  the  input  file  input.dat  in  Fig.  8  are  as  follows; 

line  1 :  Name  of  the  file  contains  input  electric  field  data  for  the  first  beam  or  circularly 

polarized  component  in  near  field  or  in  far  field. 

line  2:  Name  of  the  file  contains  input  electric  field  data  for  the  second  beam  or 

circularly  polarized  component  in  near  field  or  in  far  field. 

line  3:  Enter  ‘n’  for  near  field  data  or  ‘y’  for  far  field  data  which  is  the  output  format  of 

non2beam.F.  This  option  is  for  continuous  execution  of  the  program. 

line  4:  If  data  for  line  3  is  ‘y’,  enter  the  frame  number  to  pick  from  the  input  field  data 

which  may  contain  data  at  different  propagation  distance. 

line  5:  If  data  for  line  3  is  ‘y’,  enter  the  propagation  distance  corresponding  to  the  frame 

picked. 

line  6:  Maximum  change  in  refractive  index  due  to  nonlinearity. 

line?:  Wavelength  of  the  carrier. 

line  8:  Scale  factor  in  the  transverse  dimensions. 

line  9:  Resolution  in  x  dimension. 

line  10:  Resolution  in  y  dimension. 

line  1 1 :  Name  of  the  file  contains  the  far  field  data  for  the  first  beam  which  may  include 

data  for  all  frames  or  only  the  data  for  the  last  frame  when  output  in  HDF  images 
is  specified. 

line  12:  Name  of  the  file  contains  the  far  field  data  for  the  second  beam  which  may 

include  data  for  all  frames  or  only  the  data  for  the  last  frame  when  output  in  HDF 
images  is  specified. 

line  13:  Linear  refractive  index. 
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«  r 


line  14:  Sign  of  the  nonlinearity. 

line  15:  Length  of  the  nonlinear  medium. 

line  16:  Error  control  parameter. 

line  17:  Enter  ‘y’  for  output  data  in  HDF  format. 

line  18:  If  answer  in  line  17  is  ‘y’,  enter  name  of  the  file  containing  the  HDF  images  for 

the  first  beam. 

line  19:  If  answer  in  line  17  is  ‘y’,  enter  name  of  the  file  containing  the  HDF  images  for 

the  second  beam. 

line  20:  Number  of  output  frames  written  to  the  far  field  data  files  or  HDF  images  files. 

line  21 :  Propagation  distance  in  front  of  the  nonlinear  medium. 

line  22:  Power  ratio  of  the  first  beam  and  the  second  beam. 

line  23:  The  ratio  between  self-phase  modulation  and  cross-phase  modulation. 
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